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The investigation of spatio-temporal dynamics of bacterial cells and their molecular components requires 
automated image analysis tools to track cell shape properties and molecular component locations inside 
the cells. In the study of bacteria aging, the molecular components of interest are protein aggregates 
accumulated near bacteria boundaries. This particular location makes very ambiguous the correspon¬ 
dence between aggregates and cells, since computing accurately bacteria boundaries in phase-contrast 
time-lapse imaging is a challenging task. This paper proposes an active skeleton formulation for bac¬ 
teria modeling which provides several advantages: an easy computation of shape properties (perimeter, 
length, thickness, orientation), an improved boundary accuracy in noisy images, and a natural bacteria- 
centered coordinate system that permits the intrinsic location of molecular components inside the cell. 
Starting from an initial skeleton estimate, the medial axis of the bacterium is obtained by minimizing 
an energy function which incorporates bacteria shape constraints. Experimental results on biological 
images and comparative evaluation of the performances validate the proposed approach for modeling 
cigar-shaped bacteria like Escherichia coli. The Image-J plugin of the proposed method can be found 
online at http://fluobactracker.inrialpes.fr. 
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1. Introduction 


One of the fundamental issues addressed by Computational Cell Biology is the characterization of 


spatio-temporal dynamics of bacterial cells and their molecular components (Slepchenko, Schaff, 


et al. 2002). The rapid development of techniques for fluorescence imaging in recent years has 
opened new research opportunities, that need to cope with the availability of automated image 


analysis tools to be fully exploited. In the study of E. Coli aging (Lindner, Madden, et al. 2008), 
the molecular components of interest are protein aggregates accumulated at the old pole region of 
the aging bacterium. These subcellular components, visualized by fluorescence imaging techniques, 
need to be reliably associated to their host cells, which are generally visualized with phase-contrast 
microscopy. Figure [I] shows a pair of simultaneous frames extracted from a time-lapse phase-contrast 
and fluorescence image sequence. The automatic association between aggregates and bacteria re¬ 
quires a very accurate estimation of bacteria boundaries in this case, because aggregates accumulate 
near these boundaries (see Figure [2]). In addition to boundary accuracy, a requirement for bacteria 
modeling is the possibility to derive accurate cell measurements such as length, width, orientation, 
perimeter, etc., that are fundamental characteristics of bacteria (Osborn and Rothfield 2007). For 
instance, in bacteria like E. Coli , tracking population variability over time helps to understand 
the combinations of effects of genetic and environmental factors on cell phenotype. Finally, having 
the possibility to describe the location of subcellular components in an intrinsic cell-centered co- 
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ordinate system may be very useful to characterize intracellular dynamics (Coquel, Jacob, et al. 


2013). 



Figure 1. E. coli cells visualized through phase-contrast microscopy (left), and their molecular components (in this case 
proteine aggregates) visualized through fluorescence microscopy (right). The cell size is approximately between 3 and 5 fim. 



Figure 2. E. coli cells visualized through phase-contrast microscopy (left) and protein aggregates visualized through a fluo¬ 
rescence microscopy (right). For a better visualization, here it is shown a saturated version of the original fluorescence image. 
Protein aggregates detected by |Dimiccoli , Jacob and Moisan] ( |2015) on the original fluorescence image, are visualized as white 
points on left image. As it can be seen, they are often localized near boundaries making ambiguous the assignment to the host 
cell. The cell size is approximately between 3 and 5 fim. 


This paper aims to provide an algorithm able to extract from an image a complete parametric 
description of cigar- or rod-shaped bacteria (e.g., E. coli ) that it contains. The underlying math¬ 
ematical representation is a curvilinear skeleton (also called medial axis), which defines, after an 
appropriate (non-const ant) dilation, the boundary of the bacterium. This model is built by evolv¬ 
ing a first skeleton estimate so as to minimize an energy that promotes a good image fitting while 
enforcing the expected properties of the cell shape. In the following, we present the state of the 
art related to the use of skeletons for shape characterization and active contours, on which the 
proposed model relies. 


1.1 Background 

Segmentation, the task of partitioning an image into coherent regions, is ill-posed and simple 
continuity or homogeneity assumptions cannot cope with the large variability of object appearances. 
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In some specific applications for which one expects to see objects from a particular class, prior 
knowledge (e.g., clues on appearance or shape) may be exploited to reduce the space of variations. 
Such cues have traditionally been used by the Computer Vision community for object recognition, 


but their usefulness in segmentation has only more recently reached a common consensus (Leibe, 
Leonardis and Schiele 2006). Appearance-based methods rely on features points as Scale-Invariant 
Feature Transform (SIFT) ( Lowe|20'oT ) but their efficiency has proved to be limited since in general, 
object identity is more a function of form that a function of appearance (Siddiqi, Kimia, et al. 
2001; Siddiqi and Kimia| 1996, 1995; Leordeanu, Hebert and Sukthankar 2007). Furthermore, in 
applications such as living cell imaging, where the intensities are reduced to a minimum to avoid 
photodamage and photobleaching, key features may not be sufficiently available. 

In the use of shape as an alternative or to reinforce the use of appearance, the skeleton, defined 


by Blum (Blum 1973) as the set of medial loci of the maximal disks contained inside the object, 
has received much attention in the literature ( Bai and Latecki||2008 ; Siddiqi, Shokoufandeh, et al. 


1998; Zhu and Yuille 1995|). Compared to other shape representation such as edge fragments and 


shock graphs, the skeleton offers more flexibility in modeling spatial relationships between parts 
and has a good repeatability/distinctiveness trade-off. Although its usefulness for shape matching 
and classification on silhouettes has been recognized since the nineties, it is only recently that 


skeletons have been exploited for non-rigid object detection (Bai, Wang, et al. 2009; Trinh and 


Kimia 2011). In all these methods, the object is first recognized using reference skeletons, and then 
the skeleton is used to optimize an initial segmentation. 

The idea of optimizing an initial segmentation, under some constraints, is at the core of active 
contour (also called snakes ) methods, which specifically refine an initial approximate contour ac¬ 
cording to the image data. Early active contour model s (|Kass, Wit kin and Terzopoulos|1988 ; Osher 
and Sethian||1988 ; |Cohen [l992[ Caselles, Catte, et al.||1993 ) act by minimizing an energy function 
consisting of a data term (including image information), and a regularization term which, in gen¬ 
eral, imposes smoothness constraints on the object shape. According to the nature of the image 
information included in the data term, the existing active contour methods can be categorized into 
two types: edge-based models and region-based models. Edge-based models use local image infor¬ 
mation —typically, gradient information— ( Kass, Witkin and Terzopoulos|1988 ; Staib and Duncan 
1992; Park, Schoepflin and Kim||2001 ) to stop the contour evolution on the object boundary, while 


region-based approaches ( Cohen, Bardinet and Ayache|1992 ; Ivins and Porrill|1995 ; Zhu and Yuille 


Chesnaud, R efregi er and Bonlet||1999| |Chan and Vese 2001 ; Jeha n-Bes son, Gastand, et al. 
Zhang, Zha ng, et al.||2010 ) use global image features inside and outside the contour to con- 


1996 


2003 


trol the evolution. Region-based methods relying on local information are able to segment images 
with non-homogeneous intensities ( Li, Kao, et al.||2007 , 2008). However, as detailed by Wang et al. 
( Wang, Li, et al.||2009 ), such localization property introduces many local minima of the nonconvex 
energy functional. Consequently, the result is more dependent on the initialization of the contour. 


The initialization issue itself is since long time the object of investigations ( 

Cohen 

1991 Cohen 

and Cohen 1993;; Kass, Witkin and Terzopoulos] 1988| Bajcsy and Kovacic 1989; Jair 

i, Zhong and 

Dubuisson-Jolly||1998:; Cohen and Kimmel 1997| Amini, Weymouth and Jain 1990[ Geiger, Gupta, 

et al. 1995) and, recently, it has been partially addressed by the formulation of convex energy 


functionals (Mao, Liu and Shi 2010; Thieu, Luong, et al. 2011), for which a single global min¬ 


imum exists. According to the contour parametrization method used to model the smoothness 
constraint included in the regularization term, the existing active contour methods can be cate¬ 


gorized into three classes: level sets snakes (Osher and Sethian 1988; Caselles, Catte, et al. 1993 


Malladi, Sethian and Vemuri 1995), point-based snakes (Kass, Witkin and Terzopoulos 1988; Xu 
and Prince||1998 ) and parametric snakes ( Staib and Duncan||1992 ; Brigger, Hoeg and Unser||2000 ). 
Level sets define the 2D-curve implicitly from an evolving surface in a 3D space. They easily enable 
topology changes but are computationally expensive. The smoothness constraints are also implicit, 
that is, they are defined on the surface and not on the curve. Point-based snakes correspond to no 
parametrization or equivalently to a polygonal line (B-splines of degree zero). Parametric snakes de¬ 
fine the curve between knot points ( Brigger, Hoeg and Unser||2000[ Jacob, Bln and Unser||2004 ) via 
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a basis decomposition such as B-splines ( Marc, Menet and Medioni||1990 ; Brigger, Hoeg and Unser 
2000) or Fourier exponentials (Staib and Duncan 1992). The overall smoothness is ensured with 


curvature and eventually length minimization. When a shape-prior is introduced into point-based 


or parametric snakes, a deformable template model is obtained (Jain, Zhong and Dubuisson-Jolly 
1998). The prior shape model may be built ad hoc with analytical formulas ( Widrow||1973 ; Yuille, 
Hallinan and Cohen||1992 ; Jolly, Lakshmanan and Jain||1996 ; Lakshmanan and Grimmer [ 1996) or 
derived from a training set ( Staib and Duncan||1992||Cootes, Hill, et al.|1 993; jCootes, Taylor, et al. 


1995; Cootes and Taylor 1995; Davatzikos, Tao and Shen 2003). For analytic models, the shape 


deformations are given by the model parameters (e.g., slope for straight lines, radius for circles). 
Mostly no probability distribution is needed. Training methods, usually called active shape mod¬ 
els (Staib and Duncan 1992; Cootes and Taylor||1995 ; Davatzikos, Tao and Shen||2003 ), are more 
flexible since they can be applied easily to any reproducible shape. However, they require to define 
an average shape and deformation models, which could bias the result if, for example, the training 
set is too small. 


1.2 Active skeleton 

Instead of matching a given bacteria template to a previously computed edge map as in classical 
recognition problems, or evolving an initial contour toward the object outline as in classical active- 
contour methods, we propose an active skeleton model which evolves initial skeletal polygonal 
lines toward the true medial axes of the bacteria. The advantage of the proposed approach for 
segmentation is that it allows to introduce strong shape constraints adapted to the bacteria class 
(here, cigar-shaped bacteria), which improves the accuracy of the boundary location even in very 
noised images. 

The paper is organized as follows. In Section [2| we describe the skeletal model we are consid¬ 
ering, while Section |2.4| details the active skeleton initiation and the evolution process based on a 
variational (energy minimization) formulation. Section [3] provides results and comparisons to other 
segmentation methods. 


2. Method 

2.1 Skeleton model 

We model a digital image acquired by a contrast-phase microscopy as a real valued discrete function 
u : D C Z 2 R, where ft = [0, N] x [0, M] is a rectangle. 

Definition 1 . A ball B(x,r) centered at x and having radius r is a maximal disk of a shape J 7 if 
B{x , r) C T and VB(x' , r') C J 7 , B[x\ r') ^ B{x , r) B(x , r) (jL B(x' , r f ) 

We denote the maximal disks of a shape by the symbol B m {x , r). 

Definition 2. The morphological skeleton S G D of a bacterium B on the image u is the set of 
maximal disks centers of B : S — {x\3r > 0 : B(x,r) = B m (x,r)}. 

We represent a skeleton S by an ordered set of n points with associated radius value: S = 
{(x^, ri)} ie {i where X{ — ( Xi,yi ) G D represents the center of the maximal disk of radius r^. 
The skeleton is hence composed of n — 1 segments Sj — [xj,xjl j_i], with j G [1 ,...,n — 1]. Given 
the skeleton <S, the bacterium is built up by dilating the centers of the maximal disks according 
to their associated radius, which is linearly interpolated inside each segment. Let Sj = [xj,xf+ 1 ] a 
segment of the skeleton, then for each point x & of Sj, the corresponding radius is computed as a 
linear interpolation of rj and rj+ 1 : r&(A) = (1 — A )rj + Ar J+ i, with 0 < A < 1. 

Given an initial skeleton, the optimization of its location strongly depends on the underlying 
segment dilations and in turn to the definition of distance used for the dilation. In the following, 
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we define two distance models and, according to them, we derive the expression of the dilation 
models. 


2.2 Distance models 

In this section, we first define a simplified model of distance of a point to a segment, say d e , and 
then we define an accurate orientation based distance, say d Q , that takes the radii difference into 
account. 

Definition 3. Simplified distance of a point to a segment. Let y = (x,y) G O be any point 
of the image and s^ = [xi,xi+ 1] G S any segment of the skeleton. Denoting by p\ the projection of 
y on the line through X{ and Xi+i, p\ can be written as p\ = (1 — Ai )xi + X\Xi+\, with X\ G R. The 
simplified Euclidean-based distance of a point y to a segment Si is defined as follows. 


( \\y-Xi\\ Ai < 0 

d e (y, Si) = < \\y - x i+1 \\ Ai > 1 (2.1) 

Uly-pill A ie]o,i[ 

where || • || is the Euclidean norm. In Figj3] (up) it can be observed that the distance coincides 
with the Euclidean distance between two points only if the orthogonal projection of y on the line 
through Xi and Xi+\ lies on the segment s^. 

The value of Ai can be computed from the scalar product < y — Xi,Xi+ 1 — Xi >: \\ = 
(x-Xi)Axrf (y-yjjAyi , where Acc* = x i+ i - x i: Aj/j = y i+ 1 - y i: and Li = \J Ax* 2 + Ay; 2 (the 

length of segment sf). If the orthogonal projection p\ of y on si lies on then the distance 
of y to Si c oincides to the distance of y to p\. In this case, such a distance can be written as 
d(y,Si ) = y/(x-Xi- Ai Ax,) 2 + (y - yi - AiAyj) 2 . By introducing A e : 


( 0 Ai < 0 

A e = < 1 Ai>l (2.2) 

I Ai = (x-x^AxMy-y^Ayt Q < A X < 1 


the simplified distance from any point to a segment can be written as: d e (x,Si ) = 

\J{x - Xi - A e Axj) 2 + {y-yi- A e Ay*) 2 . 

Definition 4. Orientation-based distance of a point to a segment. 

Let y G Ll be any point of the image and Si = [xi,Xi+ 1 ] G S any segment of the skeleton. Let 
A ri = \ri+i — ri\ be the difference between the radii value associated to the extremities of the segment 
Si and Li the segment length. 

• Case 1: Li > A let ti be the common tangent to the circles with centers (xi,Xi+i) and 
radii (r^,r^ + i) respectively, lying on the same half-plane than y compared to ti. Let y be the 
orthogonal projection of y over ti, and P 2 the intersection of the line through Xi and Xi+\ 
with the line through y and y. Writing P2 as p 2 = {1 — X2 )%i + A2^+i, with X2 G R ; the 
orientation-based distance d Q is defined as follows: (see figure^ (down)): 

( || y- Xi\\ A 2 < 0 
d 0 (y,si) = < ||y-^ +1 || A 2 > 1 
[ \\y-P 2 W a 2 G]0,1[ 

where II • II is the euclidean norm. 


5 







December 23, 2016 Computer Methods in Biomechanics and Biomedical Engineering gCMBguide 


• Case 2: Li < Arv 


doil/i s i) 


\\y-Xi\\ n > n+i 
II y- x i+ ill n < n+i 


Under the assumption Li > A p 2 = (1 — A 2 )xi + \2%i+i with A 2 E R. Let A such that 
P 2 —pi = \(xi+i —Xi). Then one has simply to replace Ai with A 2 = Ai + A in the equations derived 
for the simplified distance. 



Figure 3. Euclidean distance d e (y,Si) (up) and orientation-based distance d 0 (y,Si) (down) of a point y to a skeletal segment 
Si- 


To explicit A, one can compute the angle 6 between L and Si according to sin{6 ) = 
that 


. It comes 


A =- , e A^. (2.3) 

Note again that A is an algebraic value. To deal with the case L, < Ar, one can extend the definition 
of A, for example in the following manner: 


A = 


1 - Ai A n > Li, 

-Ai - Ar* > Li, 

A =- 7 % An otherwise. 

Liy/L^-Ar? 


(2.4) 


Now one can compute in any case A 2 = Ai + A. Because of the restrictions of A 2 , the definition can 
be extended as: 


(0 A 2 < 0, 

A 0 = < 1 A 2 > 1, 

( A 2 = Ai + A A 2 g]0, 1[, 


(2.5) 
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Finally, the general expression of d Q is as follows: 

do(x, Si) = y/{x - Xi - A 0 Axj) 2 + (y- Vi- A 0 Ay*) 2 . (2.6) 


2.3 Dilation models 

Before introducing the dilations models derived from the above definitions of distance, we define 
the distance of a point to as skeleton. 

Definition 5. Distance of a point to a skeleton. Let S = {(xi,ri)}i= i v .. 5 n E D a skeleton on 
the image u : Ll —> R and y a point of D. 

d{y,S) = min |d(y,Si)} 


Definition 6. Dilation model derived from the simplified distance of a point to a seg¬ 
ment. Let S = {(xi, ^i)}ie{i..n} a skeleton and sj = [xj,Xj+ 1 ] a segment of S. For any point y Gfi, 
let Si = argmin d e (y,sj) the segment that has minimal simplified distance from y. The simplified 

j<E[l,...,n] 

dilation of S is given by the set of points such that their simplified distance to Si is inside the 
maximal disks with interpolated radius: T>(S) — {x\d e (x, Si) < (1 — A e )r^ + A e r^ + i} where A e is 
defined in equation\2.2\ 


Definition 7. Dilation model derived from orientation-based distance of a point to a 
segment. Let S — {(^, ri)} ie {i n y a skeleton and sj = [xj,Xj+ 1] the j — th segment of S. For any 
point y E Ll, let Si = argmin d 0 (y,Sj). The orientation-based dilation of S is: 




T>(S) = {y\d 0 (y. Si ) < (1 — A Q )ri + A 0 n+i} where \ Q is defined in equation 2.5. 



Figure 4. The dilation D 0 obtained by using the orientation-based distance d 0 is represented by a full red line. The dilation 
D e obtained by using the simplified distance d e differs by D 0 only in the region delimited by the two dashed vertical lines, 
where it is represented by a blue line. However, the difference between the two dilations is very small also for large difference 
radii. 


Denoting by V Q and V e the dilations obtained by using the distances d Q and d e respectively, as it 


d 0 {y,S ) =d e (y,S)fl - (^ i ) 2 -r i (A e ) 


mes, and D 0 « D e elsewhere. As it can 


can be observed on the cartoon example shown in Fig. 
in the region delimited by the vertical dashed green 
be appreciated, the dilated models differ slightly even if the difference between the radii is large. 
Since bacteria width has low variations, radii difference values are generally small compared to 
the segment length. Hence, the difference between the dilation models would be not significant. 
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Additionally, when using the simplified distance, the corresponding dilation contains less derivative 
terms and discontinuities, and consequently the variational optimization based on it would be 
computationally faster and more stable numerically. Hereafter, to simplify the notations, we will 
use the symbol d to denote either d e or and r to denote either r^(A e ) or r^(A 0 ) according to the 
desired distance model. 

2.3.0.1 General remarks. Although the distances d e and d Q are continuous, both dilation outlines 
are made up of arc circles and straight lines. Note that the scale parameter of the outline is contained 
in the radius values: big or small outlines may be described with the same skeleton. Furthermore, 
some measures of interest inherent to the bacteria class are immediate with this representation: 
orientation, thickness, perimeter and length. 


2.4 Active skeletons 


2.^.1 Skeleton initialization 

The optimization process relies on the availability of an initial skeleton for each bacterium. In 
this work, we derived it from the closed contours obtained by applying the bacteria segmentation 


method proposed in Primet, Demarez, et al. (2008). After computing the morphological skeleton 


through the use of morphological operators, a linear vectorization is performed by using an iter¬ 
ative method inspired to Wall and Danielsson (1984); Potier and Vercken (1994) but based on 


angle variations. The vectorization starts from the farthest point from the gravity center of the 
skeleton. To take into account the presence of possible holes, a spline-interpolation is used to build 
a continuous skeletal line, which is finally sampled uniformly according to the wished accuracy of 
the skeleton. 


2.4.2 Skeleton optimization 

Given an image u and the initial skeleton approximations {Sf,i = 1, ...,n}, the aim of the active 
skeleton model is to find the skeletons {S^,z = l,...,n} corresponding to the true medial axis 
of bacteria by evolving the initial skeletons. This is achieved by minimizing an energy function 
made up of two terms: a data-fidelity term and a regularity term that incorporates bacteria shape 
constraints. To handle the specific issue of several non-ovelapping objects, we propose an additional 
repulsion energy term. 


Pixel selection function 



t 

Figure 5. Pixel selection function used to express the spatial condition. For pixels inside the closed contours the signed 
difference r(x,S) — d(x,S) has an absolute value less than 1, more than 1 otherwise. Therefore, the function being symmetric, 
we have that: f ps (d — r) + / ps (r — d) = 1. 
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2.4.2.1 Data term. Images captured by a phase-contrast microscopy are characterized by a dif¬ 
ference of contrast between the foreground specimen and the background. Therefore, the data term 
has to be small when the difference of the contrast between foreground and background is large. 
Since time-lapse imaging objects contours may have large and slow variations (that is, strong low- 
frequency components), the use of a local neighborhood allows us to have a reasonable local model 
of the background. Typically, the data term corresponds to the difference between the average in¬ 
tensity inside the closed contour, and the average intensity in a ring surrounding it that corresponds 
to a high contrast between a bacterium and its local background. Since the objective function to 
be minimized has to be differentiable, the spatial condition is expressed through a smooth pixel 
selection function f ps . 


fps(t ) = 


l+sin(-|t) 

2 

0 t < -1 

1 t > 1 


*€[-1,1] 


(2.7) 


The data term E^ is expressed as the difference between the internal energy Ei n and the external 
energy E out , which are defined as follows 


Ein — 


xGl 


fp^jjXyS) ~ (j{XyS)\fc{x) 
f ps [r(x,S) - d(x,S )] 


Eout — 


E 


xGl\d(x,S)<r(x,S)+p 


fpsjdpnS) -r(x,S)]f c (x) 
f P s [d(x, S) - r(x,<S)] 


( 2 . 8 ) 


where f c is an increasing smooth contrast function allowing to adapt the imaging contrast to the 
data, p is the width of the ring surrounding the closed contours and, since x E Si C S can be written 
as x = (1 — A )xi + Aa^+i, r(x,S) — (1 — A )ri + Ar^ + i. It is easy to see that f ps [r(x,S) — d(x,S)] 
selects pixels inside the closed contour, whereas f ps [d(x,S) — r(x,S)] selects outside pixels. 


2. 4 .2.2 Regularity terms. Two regularity assumptions are embedded into the energy function by 
imposing a smooth segment angle variation and radii homogeneity. 

Most of the bacteria are weakly bent, therefore the skeleton curvature has to be as small as pos¬ 
sible. The discrete curvature at point X{ is measured by the angle of the arc at x^ that corresponds 
to the angle cq = Xf — Xi-\,Xi+\ — Xf. Numerically, we compute sm(cq) = • The 

corresponding curvature-related regularization term is: 

2 = 71—1 

E c = yi sin 2 (c>!j) (2.9) 

i =2 

Since the bacterium thickness is quite homogeneous, the difference between the radii of the 

skeleton points, {r^,i = 1, should be small. We consider the median radius value defined as 

r med — median{ri). The corresponding energy term is 
ie[l..n] 


i=n 

Eh = “ r med) 2 ( 2 - 10 ) 

2=1 

2.4- 2.3 Repulsion term. An additional repulsion term may be added when dealing with bacteria 
colonies. Denoting by d(xi,Sk) the distance between a point X{ G S\ and the skeleton <S&, the 
corresponding distance between X{ and the bacteria is t — d{xi,Sk) — r{ — rf, where r\ is the 
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interpolated radius from Sp when computing the distance d(xi,Sk)- Since t < 0 has to be avoided, 
by defining a repulsion function f rep (t ) minimal when t > 0 , the repulsion energy term E r can be 
written as: 


Er = E E frepKxi, S t ) - n - rf) (2.11) 

k Xi 

2.4-2.4 Overall energy. The global energy to minimize combines all above energy terms: E = 
aEd + bE c + cEh + dE r , where a, 6, c, d are positive weights. The energy minimization requires 
to compute derivatives according to points coordinates and radii. All derivatives are given in 
Appendix. 


3. Results 


The parameters defined in last section - a, 6, c, d, p, the contrast function f c and the repulsion 
funcion f rep - were set once for all when correctly calibrated on our dataset. The experimental 
values were: a = 10, b = 1 , c = 0.01, d — 0.1, thickness p — 2 pixels. The contrast function we 
used is f c (x) = sin 0,8 /(x), which enhances the contrast between inside and outside mean values of 
bacteria. We have chosen a quadratic repulsion function to penalize bigger overlaps: 


f (t) = / - A ) 2 if 1 < A 

Jrep\ ) S q 0 ^J ierw J ge 


(3.1) 


with A = 0.3 pixels. The algorithm was implemented in 


C with the MegaWave2 library 


3.1 Method validation 



Figure 6. The value u(x) of the pixel x is computed as a function of the distance of the pixel x from the skeletal segment Sj 
(dashed line). 


In this section, we compare our active skeleton model to the active contour model by using the 
Haussdorff distance to a synthetically generated ground truth on a set of synthetically generated 
images. Since both methods require a rough initialization we drew it by hand. To generate the 
ground truth, as well as the test data, we started by considering a bacterial colony phase-contrast 
image. First, we computed the bacteria skeletons by applying the proposed method and then, by 
relying on the skeletons obtained, we recovered a synthetic image by computing each pixel value 


1 The software in C will be made publicly available with the publication of the article 
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u(x) as 


u(x) = f 



d(x , Sj)\ 

r j( x )) 


(3.2) 


where d(x, sj ) is the distance of the pixel x to the segment Sj of a skeleton S having as extremities 
the pixels Xj and Xj+\ with associated radii Vj and r J+ i respectively (see Fig. [6]). / is the Gauss 
error function and rj{pc) is the interpolated radium computed on the projection of x on Sj {pc f in 
Fig. 6). The interpolated radium value is computed as: rj(x) = 


arj+/3rj + i 


oi+(3 


with a, /3 > 0. 



Figure 7. Synthetic image, whose noised versions constitute the test data used for validation. 


As test data, we used noised versions of the so generated synthetic image (see Figure [7]), obtained 
by adding a Gaussian white noise with increasing standard deviations. As ground truth, we used 
the boundaries from which the synthetic image has been generated. This is justified by the fact 
that, on the unnoised synthetic image, the average Haussdorff distances to ground truth of both 
methods are almost the same and of less than 1/10 of pixel (see Fig. [§]) for a standard deviation 
of noise equal to zero. 

Fig. | shows the Haussdorff distance to the ground truth of the active skeleton and the active 
contour methods. Since the active contour model assumes that the object contours are smooth, we 
previously smoothed the image by using a Gaussian kernel of standard deviation cr& = 0.5 before 
applying this method. As it is shown in Figj9j for relatively small values of the standard deviation 
of the Gaussian smoothing kernel, the performances of the active contour method slightly improves, 
since the ’’edgness” of curves is enhanced by smoothing. However, there is an upper bound to the 
level of smoothing that can be applied which is related to the scale of the image structure we are 
looking for. The proposed active skeleton model allows to enforce smoothness without becoming 
excessively sensitive to noise by integrating an a priori on the shape which is itself smooth (rod or 
cigar-shape). It is worth to remark that on real images, the image contrast across cells boundaries 
is not constant as it is for our data test images. This greatly deteriorates the performances of the 
active contour method as it can be appreciated in Fig. [lOj This Fig. shows the interest of the 
proposed method in determining the association cell-aggregate. As detailed in section [lj solving 
this association requires accurate boundary estimation since fluorescent molecular components are 
usually concentrated at the poles of bacteria. 
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Haussdorff distance to the ground truth 



Figure 8. Haussdorff distance of the active contour method (in blue) versus the Haussdorff distance of the active skeleton 
method (in red) as a function of the standard deviation of the Gaussian noise. As it can be observed, the active contour method 
is much more sensitive to the noise. 


Active contour method 



Figure 9. Haussdorff distance to the ground truth of the active contour method for increasing levels of smoothness. 


3.2 Application to the study of Escherichia Coli aging 

In Escherichia coli (E. coli) bacteria, aging-related protein aggregates accumulate at the old pole 
region of the aging bacterium. Studying the dynamics of these molecular components inside the cells 
requires automatic tools for 1) detecting protein aggregates in fluorescence images, 2) computing 
bacteria boundaries in the corresponding contrast-phase image, 3) assigning each protein aggregate 
to one cell, 4) expressing the coordinates of the protein aggregates in the basis composed of the cell 


long (the median line) and short axes (along the skeleton width). We used the method in Dimiccoli, 
Jacob and Moisan (|2015D to detect protein aggregates in fluorescence images, the method Primet, 
Demarez, et al. (2008) for segmenting the cells in contrast-phase images, and the active skeleton 
method to refine bacteria boundaries prior to the affiliation of protein aggregates to cells. Then 
the final skeletons were used in two different manners. Firstly, to localize the protein aggregates 
in the skeletal short and long axes referential. Then, the simple shape of the skeleton was used to 
estimate the total cell width and length as that of the respective skeleton. 

As a convention, we refer below to the aggregate coordinate along the long axis as the x-coordinate 
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Figure 10. (Left) Cells boundaries obtained by applying the snakes method. (Right) Cells boundaries obtained by applying 
the active skeleton method. Both used the method[Primet, Dem arez, et al.| ( |2008| > as initialization. Detected protein aggregates 
are visualized as white points. The accuracy of the active skeleton method allows a correct assignment of protein aggregates to 
cells. 



Figure 11. (a) The x— and y— coordinates of protein aggregate inside the cell correspond to the long and short axis respectively, 

(b) Histogram of the x —component of 1,644 images associated to initial trajectories. To take into account the high variability 
of cell lengths, the x-component was rescaled by division by the cell half-length, so that the cell poles are located at locations 
-1.0 and 1.0 respectively. Image adapted from the biological study [Coquel, Jacob, et al. ( |2013| ) 


and that along the short axis as the y-coordinate (see Fig. 11 (a)). In Fig. [TT] (b) is shown the 
histogram of the x —component computed over 1,644 images. It can be seen that most aggregates 
accumulate along the center, that undergoes division, and the poles of bacteria. Fig. [12] schematizes 
the increase of the cell half-length during growth that dominates the movement along the x-axis, 
also showing the time evolution of the mean displacement. 


4. Discussion 

This study has covered a number of aspects of the problem of investigating the dynamics of bac¬ 
teria cells and their molecular components. Overall the work reported produced numerical values 
for validation results and experimental results, and the discussion will be organised to highlight 
different parts of the process and to analyse the computational complexity. We will end this sec- 
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- along x-axis - along y-axis 



Figure 12. Evolution mean displacement along the x— (red) and y— axis (black). The movement along the x-axis dominated 
the cell half-length growth. Image adapted from the biological study [Coquel, Jacob, et al. ( |2013| ). 


tion by discussing the range of applicability of the proposed method. An example of initialization 
is shown on figure [l3| (a): we choose a default number of four points per bacteria (and more if 
the spline interpolation is too far from the skeleton). Even if the morphological skeleton provides 
radii values at each skeleton point, we used a constant value corresponding to the minimal contact 
radius instead (see blue outlines). In this way, the influence of the initial segmentation obtained 
by applying the method proposed by Primet el al. Primet, Demarez, et al. (2008| is minimized. 
The results correspond to what is expected: the different final skeletons are inside the bacteria 
and the implicit outlines demarcate dark areas from light outsides. Figures [T3] (b) and [l3| (c) show 
final states from initial positions given in figure [l3| (a). They highlight the interest of the repulsion 
energy E r . The second segmentation seems fully coherent with the initial image. 

From a biological point of view the segmented bacteria are too thin since they should be in 
contact according to the experimental conditions. Of course this is not the case for the above 
segmentation (e.g on figure 13 (c)) because we consider that there are bright pixels outside: the 
bright pixels isolate the bacteria from each other. For this experimental reason, the bacteria have to 
be enlarged in some way. We propose two strategies to fix this problem. The first solution consists 
in adding a positive coefficient cq n < 1 in the expression Ey = ai n Ei n — E ou t. This enables one 
to enlarge the bacteria according to the luminosity transitions as shown on figure [l3| (d). Though 
theoretically this is a good solution, the result is that different bacteria are more ore less dilated 
according to their outline intensity. Moreover, when taking low values of cq n (very big bacteria) 
the active skeleton may leak in the background. Another repulsion term from the bottom could 
be introduced but it supposes to set an additional arbitrary parameter: the minimal distance of 
bacteria to the bottom A 2 . Then the result is biased like on red dots in Fig. [l3] (e). For the 
biological study we prefer another solution. We consider that bacteria are bigger than they appear 
on the image. Then the optimization is done by considering bigger bacteria in the repulsion energy 
term E r (interaction between real bacteria) while having a thin eroded representation of bacteria 
for the data confidence energy E&. In other words, we assume that an eroded representation of the 
bacteria may fit the image. Let us call h the radius difference between real and thin bacteria. To 
adapt consequently the repulsion term E r , A in equation (3.1) is replaced by A + h to prevent 
the overlapping of the physical bacteria. The optimization process produces thin bacteria. Then to 
obtain the physical bacteria from their thin image version, a homogeneous dilation is performed 
after the optimization. In practice, it simply consists in adding h to all Ti values. Figure [l3] (f) 
presents a h = 1.5 pixels thickness difference. 

The computational complexity of the proposed method depends on three parameters: the num¬ 
ber of active skeletons, the number of skeletal points used to model each of them and the image 
resolution. More precisely, for one step of the gradient descent and for a single skeleton, the com¬ 
plexity of deriving the energy term is O(MiV), where N is the number of skeleton points and M 
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Figure 13. (a) Initial state of active skeletons (green) made up of at least 4 points with the corresponding implicit outlines 

(blue), (b) Final state of active skeletons (green) and corresponding implicit outlines (blue) without the repulsion energy E r . (c) 
Final state of active skeletons (green) and corresponding implicit outlines (blue) with the repulsion energy E r . The red skeleton 
was stopped by E r . (d) Final state of active skeletons (green) and corresponding implicit outlines (blue) with a thickness 
parameter c^ n = 0.8, still with the repulsion energy E r . (e) Final state of active skeletons with a thickness parameter a^ n = 0.6 
and an additional bottom repulsion set to A 2 = 2 pixels (effective on the red skeleton dots), (f) Final state of active skeletons 
with a repulsion of A + h = 3.3 pixels and a uniform dilation of h = 1.5 pixels. 


is the number of pixels within a box surrounding the skeleton dilation. For K skeletons without 
repulsion energy, the complexity is O(KMN). When considering the repulsion energy, it becomes 
O(KMN) + 0(iF 2 7V 2 ), where the second term corresponds to inter-bacteria distances. Since accu¬ 
rate segmentation of bacteria can be achieved with a few skeletal points, the algorithm complexity 
will mainly depend on the image resolution and on the number of bacteria in this particular ap¬ 
plication. 

We conclude by remarking that although this work was tailored to cigar- or rod-shaped bacteria, 
it would work as well for any kind of curved ovoid shape. For non-ovoid shapes, it could be adapted 
by changing the dilation model. For more complicated ramified shapes, the skeleton model itself 
should be adapted. The utility spectra of the proposed method is not restricted to bacteria, but to 
all cases in which a precise segmentation of cigar- or rod-shaped object is needed and the image 
quality is poor. For instance it could be very useful for segmenting images of insects in the study 
of insect populations, or for the segmentation and recognition of tumors and cigar- or rod-shaped 
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cells. 


Conclusions 


This paper has proposed an active skeleton approach for bacteria modeling, that presents sev¬ 
eral advantages for the study of bacteria shape properties and the dynamics of their molecular 
components. Indeed, the long and short axis of the skeleton define a reference system centered 
to the bacterium, which is crucial to study the dynamics of subcellular components inside the 
cell. In addition, the computation of shape properties such as local orientation, thickness, length 
and perimeter from the skeleton representation is straightforward. Finally, bacteria boundaries are 
computed accurately even in very noised images by introducing implicitly smooth shape priors 
into the segmentation process. The improved accuracy is important to reduce the ambiguity of the 
association molecule/cell in time-lapse fluorescence and phase-contrast imaging. 

The proposed model has been successfully used to study the dynamic of protein aggregates in 
E. Coli Coquel, Jacob, et al. (2013). Localization of proteins to specific positions inside bacte¬ 


ria is crucial to several physiological processes, including chromosome organization, chemotaxis 
or cell division and aging. The Image-J plugin of the proposed method can be found online at 
http://fluobactracker.inrialpes.fr. 

Further improvements could be discussed and developed according to the kind of data and to 
the imaging task to deal with. For example, the fidelity term of the energy functional could be 
chosen differently from the literature on data energy terms for snakes. Ad-hoc initializations could 
be derived from morphological gray-level skeletons and segment extraction. Finally, the skeleton 
model could handle more complex shapes, with more branches and angle constraints, and eventually 
with other distance definitions relaxing the circular extremity of each branch. 
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